/*Algoritmo de Borwein para calculo de PI*/

#include <stdio.h>
#include <gmp.h>
#include <time.h>
#define PRECISAO 35000000
int main(int argc, char *argv[]){
	time_t ini;
	mpf_t a, y, aux1, aux2, pi;
	int k = 0;

	mpf_init2(a, PRECISAO);
	mpf_init2(y, PRECISAO);
	mpf_init2(aux1, PRECISAO);
	mpf_init2(aux2, PRECISAO);
	mpf_init2(pi, PRECISAO);

	/*a0 = 6-4sqrt(2)*/
	mpf_sqrt_ui(aux1, 2);
	mpf_mul_ui(aux1, aux1, 4);
	mpf_ui_sub(a, 6, aux1);

	/*y0 = sqrt(2)-1*/
	mpf_sqrt_ui(aux1, 2);
	mpf_sub_ui(y, aux1, 1);

	ini = time(0);
	for(k = 0; k < 12; k++){
		/*Calculando numerador*/
		mpf_set(aux2, y);/*yk*/
		mpf_pow_ui(aux1, aux2, 4);/*yk^4*/
		mpf_ui_sub(aux1, 1, aux1); /*1-yk^4*/
		mpf_sqrt(aux1, aux1);
		mpf_sqrt(aux1, aux1);/*(1-yk^4)^(1/4)*/
		mpf_ui_sub(aux1, 1, aux1);/*1-(1-yk^4)^(1/4)*/

		/*calculando denominador*/
		mpf_set(aux2, y);/*yk*/
		mpf_pow_ui(aux2, aux2, 4);/*yk^4*/
		mpf_ui_sub(aux2, 1, aux2); /*1-yk^4*/
		mpf_sqrt(aux2, aux2);
		mpf_sqrt(aux2, aux2); /*(1-yk^4)^(1/4)*/
		mpf_add_ui(aux2, aux2, 1);


		mpf_div(y, aux1, aux2); /*y(k+1)*/
		
		mpf_add_ui(aux1, y, 1);/*(1+y(k+1))*/
		mpf_pow_ui(aux1, aux1, 4);/*(1+y(k+1))^4*/
		mpf_mul(a, aux1, a);/*ak(1+y(k+1))^4*/
		
		
		mpf_pow_ui(aux2, y, 2); /*y(k+1)^2*/
		mpf_add(aux2, aux2, y); /*y(k+1)+y(k+1)^2*/
		mpf_add_ui(aux2, aux2, 1);/*y(k+1)+y(k+1)^2+1*/
		mpf_mul(aux2, aux2, y);/*y(k+1)*(1+y(k+1)+y(k+1))^2*/
		mpf_set_ui(aux1, 2);
		mpf_pow_ui(aux1, aux1, (2*k+3));/*2^(2*k+3)*/
		mpf_mul(aux2, aux2, aux1);

		mpf_sub(a, a, aux2); /*a(k+1)*/

		mpf_ui_div(pi, 1, a);
		gmp_printf("PI(%d): %.10000000Ff\n", k, pi);
		
	}
	printf("tempo: %d segundos\n", time(0)-ini);

	mpf_clear(a);
	mpf_clear(aux1);
	mpf_clear(aux2);
	mpf_clear(pi);
	mpf_clear(y);
	
	return 0;
}	
